[a,b] = xlsread('d:\manoj\definitive_Don\lowkp_daily\daily_aver.xls');

c=unique(b);e=char(c);

d=char(b);

L = a == 99999;

a(L) = NaN;

for i = 1:105,

    L=e(i,1)==d(:,1)&e(i,2)==d(:,2)&e(i,3)==d(:,3);

    data_array = a(L,3);

    time_array = datenum(num2str(a(L,1)),'yyyymmdd');

    

    L = isnan(data_array);    

    time_array(L) = [];

    data_array(L) = [];

    

    

    b1 = min(time_array):365:max(time_array)+10;

    x = time_array'

    if length(b1) > 1,

    sp=spline(b1,data_array(:)'/spline(b1,eye(length(b1)),x(:)'));

    v=ppval(x,sp);

    

    

    plot(x,data_array);

    hold on;

    plot(b1,ones([1,length(b1)])*mean(data_array),'r.');

    plot(x,v,'g');

    title(upper(e(i,:)));

    axis([datenum(1995,1,1),datenum(2005,12,31),-inf inf]);

    datetick('x','keeplimits');

    pause;hold off;

    cla;

    end;

end;

%%



load c:\manoj\geomag\indices\kp19952006_new.mat;

load c:\manoj\projects\stormi\mat\obs_1995-2006_raw_ver_2.mat





Z_spline_removed = Z;

Y_spline_removed = Y;

X_spline_removed = X;



LL = kpdata <=100;

    

    %%

for i = 1:160,



    subplot(3,1,1);

    data_array = Z(LL,i);

    time_array = fday(LL);

    L = isnan(data_array);

    time_array(L) = [];

    data_array(L) = [];

    b1 = min(time_array):365:max(time_array)+10;

    if length(b1) > 1,

    x = time_array';

    sp=spline(b1,data_array(:)'/spline(b1,eye(length(b1)),x(:)'));

    v=ppval(x,sp);

%     plot(x,data_array);

%     hold on;

%     plot(b1,ones([1,length(b1)])*mean(data_array),'k.');

%     plot(x,v,'r');

%     axis([datenum(1995,1,1),datenum(2005,12,31),-inf inf]);

%     datetick('x','keeplimits');

%     title(['Kp <= 2' 'Z ' upper(mne(i,:))]);

%     hold off;

    end;

    

    Z_spline_removed(~L,i) = Z(~L,i) - v';

    



    

    subplot(3,1,2);

    data_array = X(LL,i);

    time_array = fday(LL);

    L = isnan(data_array);

    time_array(L) = [];

    data_array(L) = [];

    b1 = min(time_array):365:max(time_array)+10;

    if length(b1) > 1,

    x = time_array';

    sp=spline(b1,data_array(:)'/spline(b1,eye(length(b1)),x(:)'));

    v=ppval(x,sp);

%     plot(x,data_array);

%     hold on;

%     plot(b1,ones([1,length(b1)])*mean(data_array),'k.');

%     plot(x,v,'r');

%     axis([datenum(1995,1,1),datenum(2005,12,31),-inf inf]);

%     datetick('x','keeplimits');

%     title(['Kp <= 2' 'X ' upper(mne(i,:))]);

%     hold off;

    end;

    

    X_spline_removed(~L,i) = X(~L,i) - v';    

    

    subplot(3,1,3);

    data_array = Y(LL,i);

    time_array = fday(LL);

    L = isnan(data_array);

    time_array(L) = [];

    data_array(L) = [];

    b1 = min(time_array):365:max(time_array)+10;

    if length(b1) > 1,

    x = time_array';

    sp=spline(b1,data_array(:)'/spline(b1,eye(length(b1)),x(:)'));

    v=ppval(x,sp);

%     plot(x,data_array);

%     hold on;

%     plot(b1,ones([1,length(b1)])*mean(data_array),'k.');

%     plot(x,v,'r');

%     axis([datenum(1995,1,1),datenum(2005,12,31),-inf inf]);

%     datetick('x','keeplimits');

%     title(['Kp <= 2' 'Y ' upper(mne(i,:))]);

%     hold off;

    end;

    

    Y_spline_removed(~L,i) = Y(~L,i) - v';    

    

% saveas(gcf,['C:\Manoj\projects\gary\' mne(i,:)],'pdf');

display(['C:\Manoj\projects\gary\' mne(i,:) sprintf('   %d',i)]);

    %pause;

end





%%MANOJ,Your obs residual data set will be very helpful in finding out the expected contribution of external fields to magnetic navigation errors. I need that for justifying the WMMSAT requirements.

% Could you make the data available to me in a suitable ASCII format?

% Preferably:

% - One file .dat per observatory, using 3-letter code as name

% - Header: lat, lon of observatory

% - data: time (fday since 00:00:00 1/1/2000), X,Y,Z, Xres, Yres, Zres

% I need the measured field to compute the declination error

% By the way, the result can be written up as "Geomagnetic Hazard Map", so this will lead to a nice publication.

% When can I get the data?

% Stefan 



L = isnan(X);

X(L) = 99999;

L = isnan(Y);

Y(L) = 99999;

L = isnan(Z);

Z(L) = 99999;

L = isnan(X_secular_removed);

X_secular_removed(L) = 99999;

L = isnan(Y_secular_removed);

Y_secular_removed(L) = 99999;

L = isnan(Z_secular_removed);

Z_secular_removed(L) = 99999;









for i = 1:160,

filename = ['c:\manoj\projects\gary\' station_name(i,:) '.dat'];



fid = fopen(filename,'wt');



for j = 1:105192,

    fprintf(fid,'%10.4f %3.1f %10.0f %10.0f %10.0f %10.4f %10.4f %10.4f\n',fday(j),kp(j),X(j,i),Y(j,i),Z(j,i),X_secular_removed(j,i),...

    Y_secular_removed(j,i),Z_secular_removed(j,i));    

end;

fclose(fid);

fprintf('%d  %s\n',i,filename);

end;

    



%%

% Processing fo rminute data 

% 

% load c:\manoj\projects\stormi\mat\obs_1995-2006_raw_ver_2.mat X Y Z mne;

% time_array_hr = (datenum(1995,1,1): (1/(24)): datenum(2006,12,31,23,59,0));

time_array_min = (datenum(1995,1,1): (1/(24*60)): datenum(2006,12,31,23,59,0))';

S = dir('C:\Manoj\projects\gary\minute_data\*.nc');

a = sum(daysinyear([1995:2006]))*24*60;

more off;

for i = 1:length(S),



nc_fname = ['C:\Manoj\projects\gary\minute_data\' S(i).name];

fprintf('%d %s\n',i,nc_fname);   



S1 = dir(nc_fname);

if ~isempty(S1),

ncid = netcdf.open(nc_fname,'NOWRITE');

X_ID = netcdf.inqVarID(ncid,'X');

Y_ID = netcdf.inqVarID(ncid,'Y');

Z_ID = netcdf.inqVarID(ncid,'Z');

end;



x_data = netcdf.getVar(ncid, X_ID, 0, a);

x_data = double(x_data)/10;

x_data(x_data==99999.9) = NaN;



y_data = netcdf.getVar(ncid, Y_ID, 0, a);

y_data = double(y_data)/10;

y_data(y_data==99999.9) = NaN;



z_data = netcdf.getVar(ncid, Z_ID, 0, a);

z_data = double(z_data)/10;

z_data(z_data==99999.9) = NaN;

netcdf.close(ncid);



subplot(311);

plot(time_array_min,x_data);

L = isnan(x_data);

data_array = x_data(~L);

time_array = time_array_min(~L);



hold on;

b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v=ppval(time_array(1:60:end),sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v= sp(1) + sp(2) * time_array(1:60:end);

end;

%    

plot(time_array(1:60:end),v,'r');

hold off;

title(['Minute data 1995 - 2006 ' S(i).name ' X']);

datetick('x','yy');

axis([datenum(1995,1,1) datenum(2006,12,31) -inf inf]);



subplot(312);

plot(time_array_min,y_data);



L = isnan(y_data);

data_array = y_data(~L);

time_array = time_array_min(~L);



hold on;

b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v=ppval(time_array(1:60:end),sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v= sp(1) + sp(2) * time_array(1:60:end);

end;

plot(time_array(1:60:end),v,'r');

hold off;

title(['Minute data 1995 - 2006 ' S(i).name ' Y']);

datetick('x','yy');

axis([datenum(1995,1,1) datenum(2006,12,31) -inf inf]);



subplot(313);

plot(time_array_min,z_data);



L = isnan(z_data);

data_array = z_data(~L);

time_array = time_array_min(~L);



hold on;

b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v=ppval(time_array(1:60:end),sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v= sp(1) + sp(2) * time_array(1:60:end);

end;

plot(time_array(1:60:end),v,'r');

hold off;

title(['Minute data 1995 - 2006 ' S(i).name ' Z']);

datetick('x','yy');

axis([datenum(1995,1,1) datenum(2006,12,31) -inf inf]);

pause;





%saveas(gcf,['C:\Manoj\projects\gary\minute_data_plots\' S(i).name(1:3)],'pdf');

close;

end;





%%



% Processing fo rminute data 

% 

% load c:\manoj\projects\stormi\mat\obs_1995-2006_raw_ver_2.mat X Y Z mne;

% time_array_hr = (datenum(1995,1,1): (1/(24)): datenum(2006,12,31,23,59,0));

time_array_min = (datenum(1995,1,1): (1/(24*60)): datenum(2006,12,31,23,59,0))';

S = dir('C:\Manoj\projects\gary\minute_data\*.nc');

a = sum(daysinyear([1995:2006]))*24*60;

more off;

for i = 1:length(S),



nc_fname = ['C:\Manoj\projects\gary\minute_data\' S(i).name];

fprintf('%d %s\n',i,nc_fname);   



S1 = dir(nc_fname);

if ~isempty(S1),

ncid = netcdf.open(nc_fname,'NOWRITE');

X_ID = netcdf.inqVarID(ncid,'X');

Y_ID = netcdf.inqVarID(ncid,'Y');

Z_ID = netcdf.inqVarID(ncid,'Z');

end;



x_data = netcdf.getVar(ncid, X_ID, 0, a);

x_data = double(x_data)/10;

x_data(x_data==99999.9) = NaN;



y_data = netcdf.getVar(ncid, Y_ID, 0, a);

y_data = double(y_data)/10;

y_data(y_data==99999.9) = NaN;



z_data = netcdf.getVar(ncid, Z_ID, 0, a);

z_data = double(z_data)/10;

z_data(z_data==99999.9) = NaN;

netcdf.close(ncid);





L = isnan(x_data);

data_array = x_data(~L);

time_array = time_array_min(~L);



b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v_x=ppval(time_array_min,sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v_x= sp(1) + sp(2) * time_array_min;

end;

%    

L = isnan(y_data);

data_array = y_data(~L);

time_array = time_array_min(~L);





b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v_y=ppval(time_array_min,sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v_y= sp(1) + sp(2) * time_array_min;

end;



L = isnan(z_data);

data_array = z_data(~L);

time_array = time_array_min(~L);





b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v_z =ppval(time_array_min,sp);

else %data length <=1 year 

    sp=robustfit(time_array,data_array);

    v_z = sp(1) + sp(2) * time_array_min;

end;



nc_sec_fname = ['C:\Manoj\projects\gary\minute_data\secular_removed\' S(i).name(1:3) '_secular_removed.nc'];

Err = write_mag_netcdf(x_data-v_x,y_data-v_y,z_data-v_z,nc_sec_fname);



%saveas(gcf,['C:\Manoj\projects\gary\minute_data_plots\' S(i).name(1:3)],'pdf');

close;

end;





%%

% Processing fo rminute data 

% read nc files, save only sline values 1 per year X, Y Z

fid_x = fopen('C:\Manoj\projects\gary\X_spline.txt','wt');

fid_y = fopen('C:\Manoj\projects\gary\Y_spline.txt','wt');

fid_z = fopen('C:\Manoj\projects\gary\Z_spline.txt','wt');



time_array_year = [datenum(1995:2006,1,1) datenum(2006,12,31,23,59,0)]' ;

time_array_min = (datenum(1995,1,1): (1/(24*60)): datenum(2006,12,31,23,59,0))';

S = dir('D:\Manoj\gary\minute_data\orginal\*.nc');

a = sum(daysinyear([1995:2006]))*24*60;

more off;

for i = 1:length(S),



nc_fname = ['D:\Manoj\gary\minute_data\orginal\' S(i).name];

fprintf('%d %s\n',i,nc_fname);   



S1 = dir(nc_fname);

if ~isempty(S1),

ncid = netcdf.open(nc_fname,'NOWRITE');

X_ID = netcdf.inqVarID(ncid,'X');

Y_ID = netcdf.inqVarID(ncid,'Y');

Z_ID = netcdf.inqVarID(ncid,'Z');

end;



x_data = netcdf.getVar(ncid, X_ID, 0, a);

x_data = double(x_data)/10;

x_data(x_data==99999.9) = NaN;



y_data = netcdf.getVar(ncid, Y_ID, 0, a);

y_data = double(y_data)/10;

y_data(y_data==99999.9) = NaN;



z_data = netcdf.getVar(ncid, Z_ID, 0, a);

z_data = double(z_data)/10;

z_data(z_data==99999.9) = NaN;

netcdf.close(ncid);





L = isnan(x_data);

data_array = x_data(~L);

time_array = time_array_min(~L);



b1 = min(time_array):365:max(time_array)+10;

L_x = time_array_year >= min(b1) & time_array_year <= max(b1);

this_time_array_year = time_array_year(L_x);

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v_x=ppval(this_time_array_year,sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v_x= sp(1) + sp(2) * this_time_array_year;

end;

%    

L = isnan(y_data);

data_array = y_data(~L);

time_array = time_array_min(~L);





b1 = min(time_array):365:max(time_array)+10;

L_y = time_array_year >= min(b1) & time_array_year <= max(b1);

this_time_array_year = time_array_year(L_y);



if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v_y=ppval(this_time_array_year,sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v_y= sp(1) + sp(2) * this_time_array_year;

end;



L = isnan(z_data);

data_array = z_data(~L);

time_array = time_array_min(~L);





b1 = min(time_array):365:max(time_array)+10;

L_z = time_array_year >= min(b1) & time_array_year <= max(b1);

this_time_array_year = time_array_year(L_z);

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v_z =ppval(this_time_array_year,sp);

else %data length <=1 year 

    sp=robustfit(time_array,data_array);

    v_z = sp(1) + sp(2) * this_time_array_year;

end;



v_xx = zeros([1,length(time_array_year)]);

v_xx(L_x) = v_x;

v_yy = zeros([1,length(time_array_year)]);

v_yy(L_y) = v_y;

v_zz = zeros([1,length(time_array_year)]);

v_zz(L_z) = v_z;





fprintf(fid_x,'%s ', S(i).name(1:3));

fprintf(fid_x,'%8.2f ',v_xx);

fprintf(fid_x,'\n');



fprintf(fid_y,'%s ', S(i).name(1:3));

fprintf(fid_y,'%8.2f ',v_yy);

fprintf(fid_y,'\n');



fprintf(fid_z,'%s ', S(i).name(1:3));

fprintf(fid_z,'%8.2f ',v_zz);

fprintf(fid_z,'\n');





close;

end;



fclose all;











        

%% Spline fitting and plot



S = dir('E:\projects\obs_mag_data\minute\UpdateOf2011\netCDF_Absolute\*.nc');

%S = dir('E:\projects\obs_mag_data\minute\UpdateOf2011\netCDF_Absolute\new\*.nc');

a = sum(daysinyear([1995:2010]))*24*60;

more off;

%Need to work on 10

%i=74, KNZ

%for i = 1:length(S),

%for i = 42:32:74, %ESA and KNZ

for i = 55:1:55,%55   GZH

    

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30))';



nc_fname = ['E:\projects\obs_mag_data\minute\UpdateOf2011\netCDF_Absolute\' S(i).name];

fprintf('%d %s\n',i,nc_fname);   



S1 = dir(nc_fname);

if ~isempty(S1),

ncid = netcdf.open(nc_fname,'NOWRITE');

X_ID = netcdf.inqVarID(ncid,'Magnetic_Field_X');

Y_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Y');

Z_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Z');

end;



x_data = netcdf.getVar(ncid, X_ID, 0, a);

x_data = double(x_data)/10;

x_data(x_data==99999.9) = NaN;



subplot(311);



plot(time_array_min(1:60:end),x_data(1:60:end));

L = isnan(x_data);

data_array = x_data(~L);

time_array = time_array_min(~L);



clear x_data;



hold on;

b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v=ppval(time_array(1:60:end),sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v= sp(1) + sp(2) * time_array(1:60:end);

end;

%    

plot(time_array(1:60:end),v,'r');

hold off;

title(['Minute data 1995 - 2010 ' S(i).name(1:3) '.nc' ' X']);

datetick('x','yy');

axis([datenum(1995,1,1) datenum(2010,12,31) -inf inf]);



clear data_array time_array L;





subplot(312);



y_data = netcdf.getVar(ncid, Y_ID, 0, a);

y_data = double(y_data)/10;

y_data(y_data==99999.9) = NaN;



plot(time_array_min(1:60:end),y_data(1:60:end));



L = isnan(y_data);

data_array = y_data(~L);

time_array = time_array_min(~L);



clear y_data;



hold on;

b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v=ppval(time_array(1:60:end),sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v= sp(1) + sp(2) * time_array(1:60:end);

end;

plot(time_array(1:60:end),v,'r');

hold off;

title(['Minute data 1995 - 2010 ' S(i).name(1:3) '.nc' ' Y']);

datetick('x','yy');

axis([datenum(1995,1,1) datenum(2010,12,31) -inf inf]);



clear data_array time_array L;



subplot(313);



z_data = netcdf.getVar(ncid, Z_ID, 0, a);

z_data = double(z_data)/10;

z_data(z_data==99999.9) = NaN;

netcdf.close(ncid);



plot(time_array_min(1:60:end),z_data(1:60:end));



L = isnan(z_data);

data_array = z_data(~L);

time_array = time_array_min(~L);

clear z_data;

hold on;

b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v=ppval(time_array(1:60:end),sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v= sp(1) + sp(2) * time_array(1:60:end);

end;

plot(time_array(1:60:end),v,'r');

hold off;

title(['Minute data 1995 - 2010 ' S(i).name(1:3) '.nc' ' Z']);

datetick('x','yy');

axis([datenum(1995,1,1) datenum(2010,12,31) -inf inf]);

%pause;



clear  time_array_min ;

clear data_array time_array L;

saveas(gcf,['E:\projects\obs_mag_data\minute\UpdateOf2011\work\images\' S(i).name(1:3)],'pdf');

close;

end;





%% Spline fit and write the secular removed data to a netCDF file



S = dir('E:\projects\obs_mag_data\minute\UpdateOf2011\netCDF_Absolute\*.nc');



a = sum(daysinyear([1995:2010]))*24*60;



more off;



localdirectory = 'E:\projects\obs_mag_data\minute\UpdateOf2011\NetCDF_Secular_Removed\';



% Read the master station database



[master_a,master_b] = xlsread('E:\projects\obs_mag_data\minute\UpdateOf2011\work\Observatory_Master_Metadata.xls');



master_b(1,:) = [];



for i= 1:length(master_b),

master_stn(i,:) = upper(cell2str(master_b(i,1)));

end;



%define lengths of time series



length_1995_2006 = sum(daysinyear(1995:2006))*24*60;



length_1995_2010 = sum(daysinyear(1995:2010))*24*60;



where_to_store = 0;

length_of_record = a;





%for i = 1:length(S),

for i = 55:1:55, % 55 GZH



time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30))';





%read the absolute values from a netCDF file



nc_fname = ['E:\projects\obs_mag_data\minute\UpdateOf2011\netCDF_Absolute\' S(i).name];



fprintf('%d %s\n',i,nc_fname);   



S1 = dir(nc_fname);

if ~isempty(S1),

ncid = netcdf.open(nc_fname,'NOWRITE');

X_ID = netcdf.inqVarID(ncid,'Magnetic_Field_X');

Y_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Y');

Z_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Z');

end;





%create new netCDF file for writing the secular removed data



[stn , index_a ,index_b ] = intersect(master_stn,S(i).name(1:3),'rows');



metadatacellarray = master_b(index_a ,:);



lat = master_a(index_a ,1);



long = master_a(index_a ,2);



alt = master_a(index_a ,3);

    

ncid_1 = Create_Geomag_netCDF(localdirectory, lat, long, alt, metadatacellarray);



X_ID_1 = netcdf.inqVarID(ncid,'Magnetic_Field_X');

Y_ID_1 = netcdf.inqVarID(ncid,'Magnetic_Field_Y');

Z_ID_1 = netcdf.inqVarID(ncid,'Magnetic_Field_Z');



% read abs x values



x_data = netcdf.getVar(ncid, X_ID, 0, a);

x_data = double(x_data)/10;

x_data(x_data==99999.9) = NaN;



L = isnan(x_data);

data_array = x_data(~L);

time_array = time_array_min(~L);



b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp = spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v = ppval(time_array_min,sp);

else,%data length <=1 year 

    sp = robustfit(time_array,data_array);

    v= sp(1) + sp(2) * time_array_min;

end;



% write sec removed x values



x_data = x_data - v;



x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)





netcdf.putVar(ncid_1, X_ID_1, where_to_store, length_of_record,  int32 (x_data .* 10) );



clear x_data v data_array time_array L;



%read absolute y values



y_data = netcdf.getVar(ncid, Y_ID, 0, a);

y_data = double(y_data)/10;

y_data(y_data==99999.9) = NaN;





L = isnan(y_data);

data_array = y_data(~L);

time_array = time_array_min(~L);



b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v=ppval(time_array_min,sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v= sp(1) + sp(2) * time_array_min;

end;



% write sec removed y values



y_data = y_data - v;



y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)



netcdf.putVar(ncid_1, Y_ID_1, where_to_store, length_of_record,  int32 (y_data .* 10) );



clear y_data v data_array time_array L;



%read abs z vlues



z_data = netcdf.getVar(ncid, Z_ID, 0, a);

z_data = double(z_data)/10;

z_data(z_data==99999.9) = NaN;

netcdf.close(ncid);



L = isnan(z_data);

data_array = z_data(~L);

time_array = time_array_min(~L);



%fit spline



b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));

    v=ppval(time_array_min,sp);

else,%data length <=1 year 

    sp=robustfit(time_array,data_array);

    v= sp(1) + sp(2) * time_array_min;

end;



% write sec removed z values



z_data = z_data - v;



z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)





netcdf.putVar(ncid_1, Z_ID_1, where_to_store, length_of_record,  int32 (z_data .* 10) );







clear z_data v data_array time_array L;



netcdf.close(ncid_1);



clear  time_array_min X_ID X_ID_1 Y_ID Y_ID_1 Z_ID Z_ID_1 index_a index_b ncid ncid_1;







end;



%% Just plot the spline remove residuals



S = dir('E:\projects\obs_mag_data\minute\UpdateOf2011\NetCDF_Secular_Removed\*.nc');

a = (datenum(2011,1,1)- datenum(1995,1,1))*24*60;

more off;



for i = 1:length(S),



    

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30))';



nc_fname = ['E:\projects\obs_mag_data\minute\UpdateOf2011\NetCDF_Secular_Removed\' S(i).name];

fprintf('%d %s\n',i,nc_fname);   



S1 = dir(nc_fname);

if ~isempty(S1),

ncid = netcdf.open(nc_fname,'NOWRITE');

X_ID = netcdf.inqVarID(ncid,'Magnetic_Field_X');

Y_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Y');

Z_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Z');

end;



x_data = netcdf.getVar(ncid, X_ID, 0, a);

x_data = double(x_data)/10;

x_data(x_data==99999.9) = NaN;



subplot(311);



%    



plot(time_array_min(1:60:end),x_data(1:60:end) ,'b');

hold off;

title(['Minute data 1995 - 2010, SV removed  ' S(i).name(1:3) '.nc' ' X']);

datetick('x','yy');

axis([datenum(1995,1,1) datenum(2010,12,31) -inf inf]);



clear x_data;





subplot(312);



y_data = netcdf.getVar(ncid, Y_ID, 0, a);

y_data = double(y_data)/10;

y_data(y_data==99999.9) = NaN;





plot(time_array_min(1:60:end),y_data(1:60:end) ,'b');

hold off;

title(['Minute data 1995 - 2010, SV removed  ' S(i).name(1:3) '.nc' ' Y']);

datetick('x','yy');

axis([datenum(1995,1,1) datenum(2010,12,31) -inf inf]);



clear y_data;



subplot(313);



z_data = netcdf.getVar(ncid, Z_ID, 0, a);

z_data = double(z_data)/10;

z_data(z_data==99999.9) = NaN;

netcdf.close(ncid);



plot(time_array_min(1:60:end),z_data(1:60:end),'b');

hold off;

title(['Minute data 1995 - 2010, SV removed ' S(i).name(1:3) '.nc' ' Z']);

datetick('x','yy');

axis([datenum(1995,1,1) datenum(2010,12,31) -inf inf]);

%pause;



clear  time_array_min ;

clear  z_data;

saveas(gcf,['E:\projects\obs_mag_data\minute\UpdateOf2011\work\images\' S(i).name(1:3)],'pdf');

close;

end;



%% Update of 2012



% Add intermagnet definitive and preliminary. Look for definitive last

% value in the array and add preliminary after that - if any







% Definitive data



start_date = datenum(2008,1,1);

a = (datenum(2012,12,31) - datenum(2008,1,1) + 1) * 24 * 60;

time_array_min = (datenum(2008,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';

basedirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/ftp.intermagnet.org/definitive_2008_2012_nc/';



a_prelim = (datenum(2012,12,31) - datenum(2010,1,1) + 1) * 24 * 60;

time_array_min_prelim = (datenum(2010,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';

basedirectory_prelim = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/ftp.intermagnet.org/preliminary_2010_2012_nc/';



S1 = dir([basedirectory '*.nc']);

S2 = dir([basedirectory_prelim '*.nc']);



for i= 1:length(S1),

filestring(i,:) = S1(i).name(1:3);

end;



for i= 1:length(S2),

filestring_prelims(i,:) = S2(i).name(1:3);

end;



spike_stn = ['aqu';'brw';'cmo';'ded';'dmc';'eyr';'fur';'hrb';'hyb';'kiv';'lrm';'lyc';'mbo';'nur';'peg';'sit'];



[stn , index_a] = setdiff(filestring_prelims,filestring,'rows','legacy');





for i = 1:length(filestring_prelims), % all definitive data



% open preliminary data



nc_fname_prelims = [basedirectory_prelim S2(i).name];

S4 = dir(nc_fname_prelims);

if ~isempty(S4),

ncid_p = netcdf.open(nc_fname_prelims,'NOWRITE');

X_ID = netcdf.inqVarID(ncid_p,'X');

Y_ID = netcdf.inqVarID(ncid_p,'Y');

Z_ID = netcdf.inqVarID(ncid_p,'Z');

end;



x_data_prelim = netcdf.getVar(ncid_p, X_ID, 0, a_prelim);

x_data_prelim = double(x_data_prelim)/10;

x_data_prelim(x_data_prelim==99999.9) = NaN;



y_data_prelim = netcdf.getVar(ncid_p, Y_ID, 0, a_prelim);

y_data_prelim = double(y_data_prelim)/10;

y_data_prelim(y_data_prelim==99999.9) = NaN;



z_data_prelim = netcdf.getVar(ncid_p, Z_ID, 0, a_prelim);

z_data_prelim = double(z_data_prelim)/10;

z_data_prelim(z_data_prelim==99999.9) = NaN;

netcdf.close(ncid_p);







% conditional editing specific to stations



if S2(i).name(1:3) == 'ale'

    x_data_prelim(1:24*60*366) = NaN; % 2010 data is crap

    y_data_prelim(1:24*60*366) = NaN; % 2010 data is crap

    z_data_prelim(1:24*60*366) = NaN; % 2010 data is crap

elseif S2(i).name(1:3) == 'cnh'

     x_data_prelim( x_data_prelim==0 & y_data_prelim==0 & z_data_prelim==0) = NaN;

     y_data_prelim( x_data_prelim==0 & y_data_prelim==0 & z_data_prelim==0) = NaN;

     z_data_prelim( x_data_prelim==0 & y_data_prelim==0 & z_data_prelim==0) = NaN;

elseif S2(i).name(1:3) == 'mgd'

     x_data_prelim(24*60*589:end) = NaN; % > 2011.5 data is crap

     y_data_prelim(24*60*589:end) = NaN; % > 2011.5 data is crap

     z_data_prelim(24*60*589:end) = NaN; %  > 2011.5 data is crap

end;



 x_factor = 0;

 y_factor = 0;

 z_factor = 0;

 x_prelim = x_data_prelim;

 y_prelim = y_data_prelim;

 z_prelim = z_data_prelim;

 t_prelim = time_array_min_prelim;

 

 

 % spike detection and removal



x_prelim = x_data_prelim(1:end);

t_prelim = time_array_min_prelim(1:end);

if (sum(~isnan(x_prelim)) > 100 && any(intersect(spike_stn,filestring_prelims(i,:),'rows')))

    

LL = isnan(x_prelim);

x_prelim(LL) = [];

t_prelim(LL) = [];





    

[a3,b3] = robustfit(t_prelim(1:100:end),x_prelim(1:100:end));



x_prelim = x_data_prelim(1:end);

t_prelim = time_array_min_prelim(1:end);



resid = (x_prelim) - ( a3(1) + a3(2).*t_prelim);

 %LL = abs(resid) >500;

 LL = abs(resid) > 10*nanstd(resid);

 



x_prelim(LL) = NaN;

end;



%Y



y_prelim = y_data_prelim(1:end);

t_prelim = time_array_min_prelim(1:end);

if (sum(~isnan(y_prelim))  > 100 && any(intersect(spike_stn,filestring_prelims(i,:),'rows')))

LL = isnan(y_prelim);

y_prelim(LL) = [];

t_prelim(LL) = [];



[a3,b3] = robustfit(t_prelim(1:100:end),y_prelim(1:100:end));



y_prelim = y_data_prelim(1:end);

t_prelim = time_array_min_prelim(1:end);



resid = (y_prelim) - ( a3(1) + a3(2).*t_prelim);

 %LL = abs(resid) >500;

LL = abs(resid) > 10*nanstd(resid);



y_prelim(LL) = NaN;

end;

%Z





z_prelim = z_data_prelim(1:end);

t_prelim = time_array_min_prelim(1:end);

if (sum(~isnan(z_prelim))  > 100 && any(intersect(spike_stn,filestring_prelims(i,:),'rows')))

LL = isnan(z_prelim);

z_prelim(LL) = [];

t_prelim(LL) = [];



[a3,b3] = robustfit(t_prelim(1:100:end),z_prelim(1:100:end));



z_prelim = z_data_prelim(1:end);

t_prelim = time_array_min_prelim(1:end);



resid = (z_prelim) - ( a3(1) + a3(2).*t_prelim);

 %LL = abs(resid) >500;

LL = abs(resid) > 10*nanstd(resid);

z_prelim(LL) = NaN;

% 

 end;





 

[pstn,index_c ] = intersect(filestring,filestring_prelims(i,:),'rows','legacy');

 

if (~isempty(pstn))



nc_fname = [basedirectory S1(index_c).name];



S3 = dir(nc_fname);

if ~isempty(S3),

ncid = netcdf.open(nc_fname,'WRITE');

X_ID = netcdf.inqVarID(ncid,'X');

Y_ID = netcdf.inqVarID(ncid,'Y');

Z_ID = netcdf.inqVarID(ncid,'Z');

end;



x_data = netcdf.getVar(ncid, X_ID, 0, a);

x_data = double(x_data)/10;

x_data(x_data==99999.9) = NaN;



y_data = netcdf.getVar(ncid, Y_ID, 0, a);

y_data = double(y_data)/10;

y_data(y_data==99999.9) = NaN;



z_data = netcdf.getVar(ncid, Z_ID, 0, a);

z_data = double(z_data)/10;

z_data(z_data==99999.9) = NaN;



netcdf.close(ncid);





subplot(311);

plot(time_array_min,x_data);

hold on;

title(S1(index_c).name);

ylabel('X');

datetick;

subplot(312);

plot(time_array_min,y_data);

hold on;

ylabel('Y');

datetick;

subplot(313);

plot(time_array_min,z_data);

hold on;

datetick;

ylabel('Z');

 



% adjust the jumps in the preliminary data

% find the last value's index



L = isnan(x_data);

k = find(L==0,1,'last');



if (~isempty(k))

kk = find(time_array_min_prelim == time_array_min(k+1));

x_median = nanmedian(x_data(k-1000:k));

y_median = nanmedian(y_data(k-1000:k));

z_median = nanmedian(z_data(k-1000:k));

else,

    kk = 1;

    x_median = 0;

    y_median = 0;

    z_median = 0;

end;





% find the median value of the first 1000 points in the def data









% jump removal





if (~isempty(k) && length(x_prelim)> 10000 && sum(~isnan(x_prelim)) > 100 )

    

x_prelim_median = nanmedian(x_prelim(1:10000));

if isnan(x_prelim_median),

    x_prelim_median = nanmedian(x_prelim(1:100000));

end;

y_prelim_median = nanmedian(y_prelim(1:10000));

if isnan(y_prelim_median),

    y_prelim_median = nanmedian(x_prelim(1:100000));

end;

z_prelim_median = nanmedian(z_prelim(1:10000));

if isnan(z_prelim_median),

    z_prelim_median = nanmedian(x_prelim(1:100000));

end;



else

    x_prelim_median = NaN;

    y_prelim_median = NaN;

    z_prelim_median = NaN;

    

end;





if all(~isnan([x_median x_prelim_median]))

    x_factor = x_median - x_prelim_median;

else,

    x_factor = 0;

end;



if all(~isnan([y_median y_prelim_median]))

    y_factor = y_median - y_prelim_median;

else,

    y_factor = 0;

end;

if all(~isnan([z_median z_prelim_median]))

    z_factor = z_median - z_prelim_median;

else,

    z_factor = 0;

end;



end;







subplot(311);



plot(t_prelim,x_prelim + x_factor,'r');

datetick;

hold on;

title(S2(i).name);

ylabel('X');

subplot(312);



plot(t_prelim,y_prelim + y_factor,'r');

hold on;

ylabel('Y');

datetick;

subplot(313);



plot(t_prelim,z_prelim + z_factor,'r');

hold on;

datetick;

ylabel('Z');







saveas(gcf,['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/plots/' S2(i).name(1:3) '_with_spike'],'png');

display(S2(i).name(1:3));

close all;







end;





%% Validate the HDGM2013 againist a selected set of observatories


localdirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/absolute/';

S = dir([localdirectory '*.nc']);
for i= 1:length(S),
        filestring(i,:) = S(i).name;
end;

%selected_obs = ['KOU';'VSS';'TRW';'HUA'];
selected_obs = ['AAA';	'ABG';	'ABK';	'ASP';	'BDV';	'BEL';	'BFO';	...
    'BLC';	'BMT';	'BOX';	'CBB';	'CLF';	'CTA';	'CZT';	'DRV';	...
    'EBR';	'EYR';	'FCC';	'FUR';	'GNA';	'HLP';	'HRN';	'KAK';	...
    'KDU';	'KNY';	'LRM';	'LVV';	'MAW';	'MBO';	'MCQ';	'MEA';	...
    'MMB';	'NGK';	'NUR';	'NVS';	'OTT';	'PAF';	'RES';	'SOD';	...
    'SPT';	'STJ';	'THY';	'UPS';	'VIC'];



[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations

%selected observatories

[b_sel, i_aa, i_bb] = intersect(b,selected_obs,'rows');



time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a= length(time_array_min);
where_to_store = 0;
length_of_record = a;

%for i = 40:length(b_sel),
for i = 40:40,

nc_fname = [localdirectory S(i_aa(i)).name];

fprintf('%d %s\n',i,nc_fname);   

S1 = dir(nc_fname);

if ~isempty(S1),
ncid = netcdf.open(nc_fname,'NOWRITE');
X_ID = netcdf.inqVarID(ncid,'Magnetic_Field_X');
Y_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Y');
Z_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Z');
end;

% read abs x values


x_data = netcdf.getVar(ncid, X_ID, 0, a);
x_data = double(x_data)/10;
x_data(x_data==99999.9) = NaN;


L = isnan(x_data);
data_array = x_data(~L);
time_array = time_array_min(~L);

b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,

    sp = spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
    v = ppval(time_array_min,sp);
else,%data length <=1 year 
    sp = robustfit(time_array,data_array);
    v= sp(1) + sp(2) * time_array_min;

end;

subplot(311);
plot(time_array_min, v);

sp_x = sp;

clear  v data_array time_array L;

%read absolute y values


y_data = netcdf.getVar(ncid, Y_ID, 0, a);
y_data = double(y_data)/10;
y_data(y_data==99999.9) = NaN;

L = isnan(y_data);
data_array = y_data(~L);
time_array = time_array_min(~L);

b1 = min(time_array):365:max(time_array)+10;
if length(b1) > 1,
    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
    v=ppval(time_array_min,sp);
else,%data length <=1 year 
    sp=robustfit(time_array,data_array);
    v= sp(1) + sp(2) * time_array_min;
end;

subplot(312);
plot(time_array_min, v);

sp_y = sp;


clear  v data_array time_array L;

%read abs z vlues

z_data = netcdf.getVar(ncid, Z_ID, 0, a);
z_data = double(z_data)/10;
z_data(z_data==99999.9) = NaN;

L = isnan(z_data);
data_array = z_data(~L);
time_array = time_array_min(~L);
%fit spline


b1 = min(time_array):365:max(time_array)+10;
if length(b1) > 1,
    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
    v=ppval(time_array_min,sp);
else,%data length <=1 year 
    sp=robustfit(time_array,data_array);
    v= sp(1) + sp(2) * time_array_min;
end;

subplot(313);
plot(time_array_min, v);
sp_z = sp;


clear  v data_array time_array L;

eval(['save /nfs/satmag_work/mnair/projects/compare_hdgm2013/' S(i_aa(i)).name 'spline.mat' ' sp_x sp_y sp_z']);
eval(['save /nfs/satmag_work/mnair/projects/compare_hdgm2013/' S(i_aa(i)).name 'original_data.mat' ' x_data y_data z_data']);

clear x_data y-data z_data;

netcdf.close(ncid);


end;

